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ABSTRACT 

We present general relativistic hydrodynamics simulations of constant specific angular mo- 
mentum tori orbiting a Schwarzschild black hole. These tori are expected to form as a result 
of stellar gravitational collapse, binary neutron star merger or disruption, can reach very high 
rest-mass densities and behave effectively as neutron stars but with a toroidal topology (i.e. 
"toroidal neutron stars"). Our attention is here focussed on the dynamical response of these 
objects to axisymmetric perturbations. We show that, upon the introduction of perturbations, 
these systems either become unstable to the runaway instability or exhibit a regular oscillatory 
behaviour resulting in a quasi-periodic variation of the accretion rate as well as of the mass 
quadrupole. The latter, in particular, is responsible for the emission of intense gravitational ra- 
diation whose signal-to-noise ratio at the detector is comparable or larger than the typical one 
expected in stellar-core collapse, making these new sources of gravitational waves potentially 
detectable. We discuss a systematic investigation of the parameter space both in the linear and 
nonlinear regimes, providing estimates of how the gravitational radiation emitted depends on 
the mass of the torus and on the strength of the perturbation. 

Key words: accretion discs - general relativity - hydrodynamics - oscillations - gravitational 
waves 



1 INTRODUCTION 



The theory of non-geodesic, perfect fluid, relativistic tori orbiting 
a black hole has a long hist ory dating back to fundamental works 
of the late 1960s and 1970s (poyerj|l965|; |Abramowicz||l974 Fish- 
bone & Moncrief* ]HH6| r Fb*™™ w ^'' e ^i^ 1^7^ T |Kozlowskiet al 



1978). One of the most important results obtained in this series of 



investigations was the discovery that stationary barytropic config- 
urations exist in which the matter is contained within "constant- 
pressure" equipotential surfaces. Under rather generic conditions, 
these surfaces can possess a sharp cusp on the equatorial plane. 
The existence of this cusp does not depend on the choice of the 
specific angular momentum distribution and introduces important 
dynamical differences with respect to the standard model of thin 
accretion discs proposed by Shakura & Sunyaev (1973) and later 
extended to General Relativity by Novikov & Thorne (1973). 

The first important difference is that the cusp at the inner edge 
of the torus can behave as an effective Li Lagrange point of a bi- 
nary system (although this is really a circle), providing a simple 
way in which accretion can take place even in the absence of a 
shear viscosity in the fluid. The second important difference is that 
for a torus filling its outermost closed equipotential surface (the 
toroidal equivalent of a Roche lobe) the mass loss through th e cusp 
may lead to a "runaway" instability ( Abramowicz et al. 198; ). Any 



amount of material accreted onto the black hole through the cusp 
would change the black hole mass thus affecting the equipotential 
surfaces and the location of the cusp. If the cusp moved to smaller 
radial positions, the new configuration would be of equilibrium and 
no further accretion would follow. If, on the other hand, the cusp 
moved to larger radial positions, the new configuration would not 
be of equilibrium and additional material (which was previously 
in equilibrium) would accrete onto the black hole which, in turn, 
would further increase its mass. This process could trigger a run- 
away mechanism in which more and more mass is accreted onto the 
hole, evacuating the whole torus on a dynamical timescale between 
10 ms and 1 s. 

The runaway instability has also attracted attention in connec- 



tion with models of 7-ray bursts (Daigne & Mochkovitch 



1997 



Meszarosj 2002). In these models, in fact, the central engine is as 



sumed to be a torus of high density matter orbiting a stellar mass 
black hole, with intense electromagnetic emission proc esses lasting 
up to a few seconds (see Ruffert & Jankc (1995 , 2001 ) for a recent 



review on this). While many of the investigations of the runaway 
instability have concentrated on the stability properties of station- 
ary models (both in Newtonian gravity and in General Relativity), a 
time-dependent and fully relativis tic study of the runawa y stability 
has been presented only recently ( Font & Daigne 2002a). Through 
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a number of hydrodynamical simulations, Font & Daigne ( 2002a ) 



were able to show that, for constant specific angular momentum 
tori slightly overflowing their Roche lobes, the runaway instability 
does take place and for a wide range of ratios between the mass of 
the torus and that of the black hole. 

It should be noted, however, that while the instability seems 
a robust feature of the dynamics of constant specific angular mo- 
mentum tori, its existence has been severely questioned under more 
generic initial conditions. Different works, in fact, have shown that 
a more detailed modelling of the initial configurations can either 
suppress or favour the instability. Ta king into account the self - 
gravity of the torus seems to favour it ( Masuda & Eriguclv] 1997). 
The inclusion of rotation of the black hole, on the other hand, has a 



general stabilizing effect (Wilson 



1984; ^bramowicz et al 



1998) 



The same applies for tori wit h non-constant angular momen tum 
distributions, as shown first by ^Daigne & Mochkovitch (199"/) us- 
ing stationary models and by Masuda et al. (1998) with SPH time- 



dependent simula tions with a pseudo-New tonian potential. We note 
that very recently Font & Daigne ( 2002b ) have extended their rela- 



tivistic simulations to the case of non self-gravitating tori with non- 
constant angular momentum, finding that the runaway instability 
can be suppressed already with a slowly increasing specific angu- 
lar momentum distribution. A summary of the different results ob- 
tained with t he different approxima tions made so far can be found 
in Table I of |Font & Daigne| (§002a|). 

A first aim of this paper is to establish how sensitive the onset 
of the instability is on the choice of constant specific angular mo- 
mentum configurations that are initially overflowing their Roche 
lobe. To do this we adopt a mathematical and nume rical approach 
similar to the one used by Font & Daigne (2002a). It should be 



noted, however, that while we analyze the behavi our of tori that 
have masses comparable to the ones considered by Font & Daigne 



(2002a), they are also much more compact. As a result, our tori 



generically have higher rest-mass densities, in some cases almost 
reaching nuclear matter density. A second and most important aim 
of this paper is to investigate the dynamical response of these rela- 
tivistic tori to perturbations. Our interest for this has a simple jus- 
tification: because of their toroidal topology, these objects have in- 
trinsically large mass quadrupoles and if the latter are induced to 
change rapidly as a consequence of perturbations, large amounts of 
gravitational waves could be produced. 

Both of the aspects mentioned above justify in part our choice 
of terminology. Over the years, in fact, different authors have re- 
ferred to these mo dels in a number of ways, starting from the orig- 
inal suggestion of Abramowicz et al. (1978) of "toroidal stars" to 
the more recent and common denomination of "accretion tori", or 
"thick discs". Hereafter, however, we will refer to these specific ob- 
jects as tori, but also, reviving the original definition by A bramow- 
icz et al. ( 1978| ), as "toroidal neutron stars". There are three rea- 
sons for this unconventional choice. Firstly, these objects have equi- 
librium configurations with (small) finite sizes that are pressure 
supported and not accreting. In this sense, they are very different 
and should not be confused with standard accretion discs that are 
in principle infinitely extended, are generically thin because not 
pressure supported and are, of course, accreting. Secondly, these 
objects have rest-mass densities much larger than the ones usually 
associated with standard accretion discs. Thirdly, while possessing 
a toroidal topology these objects effectively behave as the more fa- 
miliar neutron stars, most notably in their response to perturbations. 

While this analogy is attractive, important differences exist be- 
tween toroidal and ordinary neutron stars, the most important being 
that toroidal neutron stars are generically unstable while spherical 



neutron stars are generically stable. More of these differences will 
appear in the following Sections. As a final remark we note that 
the idea of toroidal neutron stars might appear less bizarre when 
considering a neutron star as a fluid object whose equilibrium is 
mainly determined by the balance of gravitational forces, pressure 
gradients and centrifugal forces. In this framework then, the famil- 
iar neutron stars with spherical topology are those configurations 
in which the contributions coming from the centrifugal force are 
much smaller than the ones due to pressure gradients and gravita- 
tional forces. On the other hand, when the contributions of the pres- 
sure gradients are smaller than the ones due to the centrifugal and 
gravitational forces, a toroidal topology is inevitable and a toroidal 
neutrons^arttaij3C£CTOcs_ari^^ 
et al , ^0^2y ibr*TTecelT*s^rrmiary**ofthe*T^ 
rotating axisymmetric fluid configurations). 

The plan of the paper is as follows. In Section |^ we briefly 
recall the main properties of toroidal neutron stars, while Section H 
is devoted to a discussion of the equations to be solved and of the 
numerical approach used to solve them. After discussing our choice 
of initial data in Section^ the results of the numerical calculations 
will be presented in Section ^, reserving Section ^ to the discussion 
of the gravitational radiation produced. Finally, Section ^ contains 
our conclusions and the plans in which our future research will be 
organized. 

Throughout, we use a space-like signature ( — , +, +, +) and 
a system of geometrized units in which G = c = 1. The unit 
of length is chosen to be the gravitational radius of the black hole, 
r g = GM/c 2 , where M is the mass of the black hole. When useful, 
however, cgs units have been reported for clarity. Greek indices are 
taken to run from to 3 and Latin indices from 1 to 3. 



2 ANALYTIC SOLUTIONS FOR STATIONARY 
CONFIGURATIONS 

In what follows we recall the basic properties of stationary toroidal 
fluid configurations in a curved spacet ime and the interested r eader 
will find a more detailed discussion in Font & Daigne (2002a). The 
considerations made here will be useful only for the construction of 
the background initial model which we will then perturb as detailed 
in Section 0. 

Consider a perfect fluid with four- velocity u and described by 
the stress-energy tensor 



rpfl 



(e + p^u^u 1 + pg^ v = phii^u" + pg p 



(1) 



where g^ v are the coefficients of the metric which we choose to 
be those of a Schwarzschild black hole in spherical coordinates 
(t, r, 8, (f>). Here, e, p, p, and h = (e + p)/p are the proper en- 
ergy density, the isotropic pressure, the rest mass density, and the 
specific enthalpy, respectively. In the following we will model the 
fluid as ideal with a polytropic equation of state (EOS) p = np~' = 
pe(7 — 1), where e = e/p — 1 is the specific internal energy, k 
is the polytropic constant and 7 is the adiabatic index. Also, for 
simplicity we will consider the fluid not be magnetized. This may 
represent a crude approximation given that toroidal neutron star are 
probably created by material originally magnetized and that very 
large magnetic fields can be easily produced when rapid sh earing 
motions are present in high ly conducting magnetized fluids ( Spruit 



1999; Rezzolla et al 



2000| ). 



The fluid is assumed to be in circular non-geodesic mo- 
tion with four-velocity 11° = (it', 0, 0, u*) = u*(l,0,0,n), 
where Q = fi(r,0) = is the coordinate angular velocity 
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as observed from infinity. If we indicate with I the specific an- 
gular momentum^ (i.e. the angular momentum per unit energy) 
t = —u^/ut, the orbital velocity can then be written in terms 
of the angular momentum and of the metric functions only, as 

n = -t(.9tt/9<i»t>)' 

From the normalization condition for the four- velocity vector, 

u a u a = —1, we derive both the total specific energy of the fluid 

element, u t and the redshift factor, it* as 



M 2 = 



g<t><t>9tt 



ut = -- 



gtt + fl 2 g^ 



(2) 



Under these assumptions, the equations of motion for the 
fluid can be generically written as /i l ^V M T M ^ = where h^ v = 
9p,v + u^u v is the projector tensor orthogonal to u and V the co- 
variant derivative in the Schwarzschild spacetime. Enforcing the 
conditions of hydrostatic equilibrium and of axisymmetry simpli- 
fies the above equations considerably. Furthermore, if the contribu- 
tions coming from the self-gravity of the torus can be neglected, 
the relativistic hydrodynamics equations reduce to Bernoulli-type 
equations 



e + p 



i - m 



(3) 



where i = r,9 and W — W(r,8) = ln(wA 

The simplest solution to equations (§) is the one with £ = 
const., since in this case the equipotential surfaces can be com- 
puted directly through the metric coefficients and the value of the 
specific angular momentum. Note that at any point in the (r, 9) 
plane, the potential W can either be positive (indicating equipo- 
tential surfaces that are open) or negative (indicating equipoten- 
tial surfaces that are closed). The case W = refers to that spe- 
cial equipotential surface which is closed at infinity. Interestingly, 
closed equipotential surfaces contain local extrema and in the equa- 
torial plane these mark two very important points. There, in fact, 
ViW = = Vip and an orbiting fluid element would not experi- 
ence any net acceleration, with the centrifugal force balancing the 
gravitational one exactly. These points correspond to the (radial) 



positions of the cusp, 



. and of the "centre" of the torus, 



. At 



these radial positions the specific angular momentum must be that 
of a Keplerian geodesic circular orbit 



I 2 = f K 



Mr 



(r-2M) 2 



(4) 



which can effectively be used to calculate the position of both the 
centre and the cusp. 

In the case of a torus with constant specific angular momen- 
tum, it is straightforward to show that the position of the cusp is 
located between the marginally bound circular orbit, rt,, and the 
marginally stable circular orbit, r ms , of a p oint-like particle orbit- 
ing the black hole ( Abramowicz et al. 197S ). Note that the position 



of the inner edge of the torus n n and the position of the cusp r cusp 
need not coincide and indeed r; n can be chosen to be anywhere be- 
tween the cusp and the centre. Once made, however, the choice for 
rin also determines the position of the outer edge of the torus on the 



1 Note that this is not the only definition for the specific agular momentum 
used in the literature. Often, in fact, the specific angular momentum is de- 
fined as £' = Uj, because this is a constant of geodesic (i.e. zero-pressure) 
motion in axially symmetric spacetimes. When the pressure is non-zero, 
on the other hand, hu^ is a constant of motion, while £' is not. For axi- 
ally symmetric, stationary spacetimes £ = —u^/ut is constant for both 
geodesic and perfect fluid motion. 



equatorial plane through the constraint that both points belong to 
the same equipotential surface, i.e. W(r ou t, it/ 2) — W(r- m , tt/2). 

A particularly attractive feature of tori with constant specific 
angular momentum is that once the background spacetime and the 
value of the specific angular momentum have been fixed, the an- 
gular velocity Q. = Q(r,6) is fully determined. Furthermore, if 
a polytropic EOS is used, the Bernoulli equations ^ can be in- 
tegrated analytically to yield the rest-mass density (and pressure) 
distribution inside the torus as 



p(r,9) 



7-1 

KJ 



[exp(Wm - W) - 1] 



1/(7-1) 



(5) 



where Win = W(ri n , tt/2). Once the rest-mass distribution is 
known, the total rest-mass of the torus can be easily calculated as 



M t , 



/ t j3 

psj —gu a x , 



(6) 



where yf— g = r 2 sin 9, while the total mass-energy in the toroidal 
neutron star is computed as 



M t 



(t; + tI + t! - t*) 



-gd 3 x . 



(7) 



where d s x = drd9d(j> is the coordinate volume element. (Note 
that for simplicity, hereafter, we will refer to the mass-energy of 
the toroidal neutron star as the "mass" of the toroidal neutron star.) 

Depending on the value of £ chosen and in particular on how 
this compares with the specific angular momenta corresponding to 
orbits that are marginally bound, I m b = 4, or that are marginally 
stable, £m B = 3\/3/2, different configurations can be built. A de- 
t ailed classification of the models can be found in the literature 
( Abramowicz et al]|l978| ; Font & Daigne 2002a); here we simply 
recall that if £ ms < I < £ m b, there will be an equipotential sur- 
face closed at a finite radius and possessing a cusp. As a result, a 
stationary toroidal neutron star of finite extent can be built and this 
will represent our fiducial unperturbed toroidal neutron star. 



3 MATHEMATICAL FRAMEWORK 
3.1 Hydrodynamic equations 

To preserve the conservative nature of the equations of general rel- 
ativistic hydrodynamics, namely the local conservation of baryon 
number and energy-momentum, it is convenient to cast them in the 
form of a flux-conservative hyperbolic system through the intro- 
duction of suitable "conserved" variables rather than in terms of 
the ordinary fluid, or " primitive", variables. In this case, the equa - 



tions assume the form (Banyuls et al 
<9U(w) d[aF r {w)] d[aF e {w)] 



1997; Font&Ibanez 



1998a) 



dt 



+ 



dr 



+ ■ 



09 



S(w) 



(8) 



where a = y/— goo is the lapse function of the Schwarzschild met- 
ric and where U(w) = (D, S r , So, S<f,) is the state-vector of the 
evolved variables. The other vectors F s and S appearing in ^ rep- 
resent the fluxes and sources of the evolved quantities, respectively. 
The relation between the conserved and primitive variables in the 
vector w = (p, Vi, e) are given through the following set of equa- 
tions 



D 



phT 2 Vj 



(9) 



supplemented with the ideal-fluid EOS. Note that the covariant 
components of the three- velocity are defined in terms of the spatial 
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3-metric jij to be Vi 



where v 1 = u 1 /au l . (Although in 



axisymmetry, we evolve also the azimuthal component of the equa- 
tions of motion, so that the index j in Eqs. (^) takes the values j = 
r, 9, 4> ) The Lorentz factor V measured by a local static observer 
and appearing in Eqs. (^) is defined as V = au l = (1 — v' 2 )~ 1 ^ 2 , 
with v 2 = 'yijv'v-' . The specific expressions for the components of 
the flux v ectors F' and of the source vector S can be found in F ont 
& Daigne" p002^ )! 



3.2 Spacetime evolution 

The general relativistic hydrodynamics equations we solve as- 
sume that the fluid moves in a curved spacetime (provided by the 
Schwarzschild solution) that is static. The onset and development 
of the runaway instability, on the other hand, depends crucially on 
the response of the fluid to variations of the spacetime and in par- 
ticular of its longitudinal part. To follow this in a self-consistent 
manner would require the solution of the Einstein field equations 
together with those of relativistic hydrodynamics. Computationally, 
however, this is a much harder problem. Fully numerical relativity 
codes evolving black hole spacetimes with perfect fluid matter (ei- 
ther in two or three spatial dimensions) are being developed only re - 



cently ( Brandt et al 



200C 



Shibata & Uryu 2000; Pont et al. 2002) 



In addition to this, the need of high grid resolution and of computa- 
tional timescales that are much larger than the dynamical one, may 
still be too taxing for full numerical relativity codes. 

To avoid the solution of the full Einstein equations and yet 
simulate the onset and d evelopment of t he inst ability we follow 
the approach proposed by Font & Daigne (2002a). Most notably, at 



each timestep we calculate the accretion rate at the innermost radial 
point r m i n of the grid as 



m(r min ) = -2tt 



-gDv r d9\ 



(10) 



and thus determine the amount of matter accreted onto the black 
hole as 



M n+1 =M * +At7h n ( rinlll ) 



(ID 



where the upper indices refer to a given time-level. Once the new 
mass of the black hole has been computed, the relevant metric func- 
tions are instantaneously updated as 



(12) 



so that p M „ describes the spacetime at the time-level n + 1, over 
which the hydrodynamical equations will be solved. We note that 
to be consistent the transfer of angular momentum from the torus 
to the black hole should also be taken into account. While we have 
not considered this here, the interested reader will find the details 
on a procedure t o account for the angula r momentum transfer onto 
the black hole in iFont & Daignel (Eo02bh. 



Our approach for the spacetime evolution is clearly an approx- 
imation and it masks important features such as the response of the 
black hole to the accreted mass and the corresponding emission 
of gravitational radiation. Nevertheless, this approximation is of- 
ten very good especially when the toroidal neutron stars are not 
very massive and the rest-mass accretion rates are therefore small. 
In these cases, then, the fractional variation of the black hole mass 
between two adjacent time-levels is minute and the spacetime evo- 
lution can be treated effectively as a discrete sequence of stationary 
spacetimes. 



3.3 Numerical approach 

The numerical code used in our comp utations is bas e d on a 
code that as has been firs t descri bed in Font & Ibanez (1998a) 
(see also Font & Daigne (2002ah). This code performs the nu- 
merical integration of system d8|) using upwind high resolution 
shock-capturing ( HRSC ) schem es based on approximate Riemann 
solvers (see, e.g. Fonl (200() and references therein). Exploit- 
ing the flux conservative form of equations (^), the time evo- 
lution of the discretized data from a time-level n to the subse- 
quent one n + 1 is performed according to the following scheme 



At 
Ar 

At 
A8 

+ At Si 



•l/2,j 



i,3+l/2 



(13) 



where the subscripts i,j refer to spatial (r,9) grid points, so that 
U(r,:, 9j t"). The inter-cell numerical fluxes, F' i±lj/2 ■ 

are computed using Marquina's approximate Rie- 



i,j±l/2- 



U" 

and F 

mann solver (Donat & Marquina, 1996). A piecewise-linear cell 
reconstruction procedure provides second-order accuracy in space, 
while the same order in time is obtained with a conservative two- 
step second-order Runge-Kutta scheme applied to the above time 
update. 

Our computational grid consists of iV r x Ne zones in the ra- 
dial and angular directions, respectively, covering a computational 
domain extending from r m i n = 2.1 to r max = 30 and from to 
■k. We have used numerical grids with different number of zones, 
finding that the truncation error is reduced to satisfactory values 
when 7V r = 250 and Ne — 84. All of the results presented in the 
paper have been computed with this number of grid-points. The ra- 
dial grid is logarithmically spaced in term of a tortoise coordinate 
r* = r + 2M ln(r /2M — 1), with the maximum r adial resolution 
at the innermost grid being Ar = 6 x 10~ 4 . As in Font & Daigne 



(2002a), we use a finer angular grid in the regions that are usually 



within the torus and a much coarser one outside. A grid-point be- 
longs to the external surface of the initial unperturbed torus when 
Ut(r, 9) = u t , in = it* (fin, 7r /2). This equation defines the merid- 
ional section of the surface as a closed polar curve C of equation 



sin 9 - 



<i/(l-2M/r) 



r2 ( U lin 



■ 2M/r) 



1/2 



(14) 



The angular extension 9 m of the unperturbed torus can then be 
computed by searching for the local extrema of the curve C. As 
a result, in most of our simulations 75% of the angular grid points 
are uniformly distributed in the range [9 m , t — 9m], while the re- 
maining points cover the external region. 

The boundary conditions adopted, the treatment of the vac- 
uum region outside the torus with a low density atmosphere, and 
the procedure for recovering physical variables fro m the conserved 
quantiti es D and Si are the same as those used by Font & Daigne 
(2002a). The interested reader is referred to that work for further 



details. 
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4 INITIAL DATA 

Simulating a dynamical instability with a numerical code brings up 
the problem of suitable initial conditions. A natural choice would 
be that of a configuration that is in equilibrium, where the latter 
is just a marginally stable one. In this case, then, any perturbation 
would move the configuration away from the equilibrium, inducing 
the instability on a finite timescale. While we see this happen reg- 
ularly in Nature, it is rather difficult to simulate it numerically, the 
major obstacle being the need of performing the numerical simula- 
tions on those (short) timescales that can be afforded computation- 
ally. Fortunately, however, there are ways of by-passing this limita- 
tion and these generally consist of choosing an initial configuration 
which is already slightly out of equilibrium. By controlling the de- 
viation away from the equilibrium in some parametrized form, the 
timescale for the development of the instability can then be reduced 
to values that are compatible with the computational ones. 

An approach of this type has been used in the past also to sim- 
ulate the runaway instability, where a measure of deviation away 
from the unstable equilibrium was made in terms of the potential 
difference A Win > at the inner edge of the torus. This quantity, 
defined as A Win = Win — W CU sp, accounts for the potential jump 
on the eq uatorial plane between the inner edg e of the torus and 
the cusp ( Igumenshchev & Beloborodov 199'/). By simply vary- 
ing the value of AWi n , it is then possible to select a configuration 
corresponding to a torus inside its Roche lobe and for which no 
mass transfer is possible (i.e. AWi n < 0), or a torus overflow- 
ing its Roche lobe and therefore accreting onto the black hole (i.e. 
AWin > 0). The case limiting the two classes of solutions, (i.e. 
AWin = 0) refers to a configuration that is just marginally sta- 
ble to the runaway instability, which will therefore develop over 
an infinite timescale. (Note that this condition is also equivalent to 
setting r in = r cusp .) 



All of the models considered by Font & Daigne ( 2002a ) have 



been constructed with potential differences AWin > 0, so that the 
outermost potential surface is not closed at the cusp but reaches the 
black hole. After truncating the torus at r = r CUS p, the simulations 
were carried out by evolving the set of equations discussed in Sec- 
tion 3.1. With this choice, a small fraction of the initial fluid con- 



figuration (i.e. all the one residing outside the Roche lobe) is out of 
equilibrium. Of course, this is n ot the only way of tr i ggerin g the in- 
stability. In the calculations by Vlasuda & Eriguchj (19971 the size 
of the torus was expanded by a small amount to overflow its Roche 
lobe and to set the configuration out of equilibrium. While it has 
been argued that the occurrence of the instability is not much af - 
fected by the choice of the initial model (Masuda & Eriguchi 1997) 
we have here followed a different approach to the problem of initial 
condition for the runaway instability as discussed below. 



4.1 Introducing a perturbation 

An important difference with respect to earlier works in our pre- 
scription of the initial data is that we have considered models with 
a potential barrier AWi n < 0. As a result, these represent con- 
figurations that are either marginally stable (i.e. AWin = 0) or 
even stable (i.e. AWi n < 0) with respect to the runaway insta- 
bility. Since these configurations cannot develop the instability on 
a finite timescale, we have introduced parametrized perturbations 
that would induce a small outflow through the cusp. More specif- 
ically, we have modified the stationary equilibrium configuration 
discussed in Section ^ with a small radial velocity which we have 
expressed in terms of the radial inflow velocity characterising a rel- 



ativistic spherically symmetric accretio n flow onto a S chwarzschild 
black hole, i.e. the Michel solution ( Michej 1972). Using r\ to 
parametrize the strength of the perturbation, we have specified the 
initial radial (covariant) component of the three- velocity as 



v r = 



(15) 



We regard this choice of initial data as a more realistic one for 
at least two reasons. Firstly, in this way only a small region of the 
fluid configuration, (i.e. the one located near the cusp) is effectively 
out of equilibrium since v r falls for large radii. Secondly and more 
important, an initial configuration of this type is much closer to the 
one that might be produced in Nature. We recall, in fact, that tori of 
the type considered here are expected to form in a number of differ- 
ent even ts such as the collapse of supcrmassive neutron stars (V ietri 



& Stella 1998), or th£iron^rare^olka£se_of^^ 
Fayden & Woosley' TIl)9*9] )"'o7h*e71lcenar^^ 

objects involve the coalescence of a binary system, either consist- 
ing of two neutron stars (especially if they have unequal masses, 
Shibata 2002) or consisting of a black hole and of a neutron star 
which is then disru pted by the intense tidal field ( Lee & Kluzniak 



1999a b 



Lee 



200C). In all of these catastrophic events, the newly 



formed torus will be initially highly perturbed and is expected to 
maintain also a radial velocity in addition to the orbital one . In re - 
cent Newtonian simulations performed by Ruffert & Janka (1999) 



the torus resulting from the dynamical merging of two neutron stars 
was observed to oscillate and accrete onto the newly formed black 
hole. The average inflow velocity in the central region of the newly 
formed torus was measured to be ~ 3 x 10 -3 , whereas at very 
small distances from the black hole the fluid was infalling much 
more rapidl y. To be consistent with the estimates provided by R uf- 
fert & Janka ( 1999[ ), we have chosen the parameter rj in the range 
[0.001, 0.06], corresponding to an average inward radial velocity in 
the range [0.0002,0.04], respectively. However, simulations with 
values as small as r\ = and as large as r\ — 0.17 have also 
been performed, but these have not introduced qualitatively new 
features. It should also be noted that because the orbital velocities 
are at least one order of magnitude larger than the radial ones in- 
duced through the perturbations, the contribution of the latter to the 
kinetic energy budget is rather small even when large values of rj 
are considered. 

An aspect of our initial models worth underlining is that while 
in principle the mass flux should be zero when AWin = 0, this 
is never the case numerically. The unavoidable mass flux induced 
at the cusp can be made arbitrarily small after a suitable choice of 
the mass of the torus and of the strength of perturbation. This rep- 
resents an important possibility because in the case of very small 
rest-mass accretion rates, the variations in the spacetime metric can 
be neglected and we can therefore investigate the response of the 
toroidal neutron stars in the absence of metric variations. We will 
refer to this regime as the one with a "fixed" spacetime to dis- 
tinguish it from the "dynamical" spacetime regime, in which the 
accretion rate is not negligible and metric functions need to be up- 
dated following the procedure of Section 3.2. 



With the choice of initial conditions discussed above, we have 
evolved a large number of models covering only part of the rel- 
evant parameter space. The properties of the different models are 
summarized in Table which contains the ratio between the mass 
of the toroidal star and of the black hole, the polytropic constant, 
the specific angular momentum, as well as all the relevant radii of 
the tori. Each of the models in Table |l| has been simulated for at 
least four different values of the parameter r), both on a fixed and 
on a dynamical spacetime. In all of the simulations we have kept 
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Table 1. Main properties of the constant angular momentum toroidal neutron star models used in the numerical calculations. From left to right the columns 
report: the name of the model, the star-to-hole mass ratio M&/M, the polytropic constant n, the specific angular momentum £ (normalized to M), the inner 
and outer radii of the toroidal neutron star r; n and r ou t, the radial position of cusp r cusp , the radial position of the centre r cen t re (all radii are in units of the 
gravitational radius r g ), and the orbital period at the centre of the torus t or b, expressed in milliseconds. The last two columns indicate the density at the centre 
of the torus and the average density of each model, respectively, both in cgs units. All of the models share the same mass for the black hole, M = 2.5Mq and 
adiabatic index 7 = 4/3. 



Model 


Mt/M 


K 


i 


'in 


r ou t 


'"cusp 


^centre 


^orb 


pecntre 


<P> 






(cgs) 












(ms) 


(cgs) 


(cgs) 


(a) 


1. 


4.46xl0 13 


3.8000 


4.576 


15.889 


4.576 


8.352 


1.86 


1.14X10 14 


4.72X10 12 


(b) 


0.5 


5.62X10 13 


3.8000 


4.576 


15.889 


4.576 


8.352 


1.86 


5.72X10 13 


2.36xl0 12 


(c) 


0.1 


0.96xl0 14 


3.8000 


4.576 


15.889 


4.576 


8.352 


1.86 


1.14X10 13 


4.73X10 11 


(d) 


0.05 


1.21X10 14 


3.8000 


4.576 


15.889 


4.576 


8.352 


1.86 


5.73 xlO 12 


2.36X10 11 


(e) 


0.1 


7.0xl0 13 


3.7845 


4.646 


14.367 


4.646 


8.165 


1.80 


1.61 xlO 13 


6.43 xlO 11 


(f) 


0.1 


l.OxlO 14 


3.8022 


4.566 


16.122 


4.566 


8.378 


1.87 


LlOxlO 13 


4.48X10 11 


(g) 


0.1 


2.0xl0 14 


3.8425 


4.410 


21.472 


4.410 


8.839 


2.03 


4.96xl0 12 


2.01 xlO 11 


(h) 


0.1 


3.5 xlO 14 


3.8800 


4.290 


29.539 


4.290 


9.246 


2.17 


2.41 xlO 12 


8.12X10 10 



fixed: the adiabatic index (taken to be that of a degenerate relativis- 
tic electron gas) 7 = 4/3, and the initial black hole mass, which 
we have chosen to be M = 2.5Mq for comparison with the results 
of^ont&Daigne| (|2002a|). 



5 NUMERICAL RESULTS 

In what follows we will discuss in detail the dynamics of the per- 
turbed toroidal neutron stars summarized in Table |l|. In particular, 
we will first report the results about the runaway instability and 
subsequently will discuss the long-term dynamics of the toroidal 
neutron stars in response to perturbations. 



5.1 Dynamical spacetime: the runaway instability 

Since the code used in our simulations is similar but distinct from 



the one employed by Font & Daigne ( 2002a ), we have first tested it 
against the results published by these authors. More precisely, we 
have considered out-of-equilibrium initial conditions (AWi n > 0), 
evolving this configuration using the set of equations (^) and the 
metric update (|l2|). The results obt ained agree (wi th differences 
less than one percent) with those by Font & Daigne (2002a), indi- 
cating that with this choice of initial conditions the runaway insta- 
bility develops rapidly, on timescales that are progressively smaller 
as the mass ratio Mt/M and the initial potential jump, A Win, 
are increased. The occurrence of the instability is signalled by the 
exponential growth of the rest-mass accretion rate which rapidly 
reaches super-Eddington values (cf. Figure 8 of Font & Daigne, 
2002a). It should be noted that while the simulated accretion rates 
are many orders of magnitude larger than the Eddington limit (this 
is ~ 1.2 x 10 -16 Mq/s for the black hole considered here), these 
mass fluxes are also the ones required to account for the large ener- 
getic release observed in 7-ray bursts. 

After this validating test we have investigated the onset and 
development of the runaway instability using the initial conditions 
discussed in Section ^| Shown in Figure [l] is the evolution of the 
rest-mass accretion rate for model (a) (see Table |j|) and with three 
different values of initial velocity perturbation, 77. The time is ex- 
pressed in terms of the orbital period t orD = 2-7r/f2 C entre of the 
centre of the toroidal neutron star and is reported in Table [l] for 
the different models considered. Note that the minimum rest-mass 
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Figure 1. Evolution of the rest-mass accretion rate for model (a) and for 
different values of 77 when the spacetime is allowed to vary. The time is 
expressed in terms of the timescale for the centre of the toroidal neutron 
star to perform an orbit. The inset focuses on the case in which 77 = 0.06, 
showing the results on a logarithmic scale. 



accretion rate in Figure |l| is never zero but ~ O.O2M0 /s, even 
initially. This is just the cumulative result of the tenuous ambient 
atmosphere, which is always producing a tiny and constant in time 
mass overflow at the cusp, coupled with the use of very high density 
matter which amplifies the effect. 

The behaviour of the mass flux reported in Figure [l| incor- 
porates two important features. Firstly, it shows that the runaway 
instability does occur also with this choice of initial data and that 
the growth-rate is shorter for larger initial velocity perturbations 
(i.e. for larger values of 77). As mentioned in the Introduction, this 
is an important point confirming that Roche lobe overflowing is not 
a necessary condition for the development of the instability, at least 
in constant specific angular momentum tori whose self-gravity is 
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Figure 2. Rest-mass accretion rate for two different mass ratios in a dynam- 
ical spacetime evolution. The solid and dotted lines correspond to models 
(b) and (c), respectively, and a perturbation with r) = 0.06 was used. 

not considered. The inset of Figure [l] shows the simulation with 
77 = 0.06 on a logarithmic time scale, and should be considered 
the evolution that more than any other of our sampl e resembles the 
one ob served for Roche lobe overflowing tori by Font & Daigne 
(b002ab, 



The behaviour observed in Figure [j] with different values of 
the perturbation velocity has a simple interpretation. In order to 
trigger the instability, in fact, a certain fractional change in the 
mass of the black hole (and therefore in the spacetime curvature) 
needs to be reached. If the initial perturbation is large, a consider- 
able amount of matter is accreted onto the black hole already after 
the first oscillation in the accretion rate and the instability is there- 
fore able to develop very rapidly. If, on the other hand, the strength 
of the perturbation is small, much less mass will be accreted dur- 
ing each oscillation and many more will be needed to produce the 
fractional change in the black hole mass that will accelerate the de- 
velopment of the instability. In this case, the instability will develop 
on larger timescales, which do not scale linearly with r\. 

The second novel feature to notice in Figure |l| is that, after 
the system has relaxed from the initial conditions (i.e. at about 
t/t OI b — 2) the secular growth in the rest-mass accretion rate is ac- 
companied by an oscillatory behaviour with increasing amplitude. 
These oscillations are less regular in the case of the high ampli- 
tude perturbation (i.e. for r\ = 0.06) but are much more regular 
as the strength of the perturbation is gradually reduced. It is im- 
portant to notice that when very low amplitude perturbations are 
used, the amount of accreted matter is so small that new features 
can be revealed. In the case of the simulation for r\ — 0.02 in 
Figure [l] (dashed line), for instance, it is possible to distinguish 
at least three different stages. Most notably, an initial stage for 
t/t OT b < 18' during which the rest-mass accretion rate is very 
small and does not manifest a regular oscillatory behaviour. De- 
spite the apparent quiescence, and as it will become clearer in the 
following Section, during this stage the toroidal neutron star is not 
at all static and other hydrodynamical quantities manifest a differ- 
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Figure 3. Evolution of the rest-mass of the toroidal neutron star (normal- 
ized to its initial value) in a dynamical (continuous line) and in a fixed 
(dashed line) spacetime simulation. The calculations refer to model (a) with 
an initial perturbation r\ = 0.06. 



ent behaviour. This stage is then followed by an intermediate stage 
for 18 < t/t OI h < 45, during which the rest-mass accretion rate 
shows the secular growth already observed for higher amplitudes 
perturbations. While this happens, the toroidal neutron star period- 
ically enters short phases during which essentially no accretion is 
present (i.e. when the accretion rate reaches its minimum). Even- 
tually, a third final stage sets in for t/t OI b > 45, during which the 
instability starts to develop more clearly. This can be deduced from 
the fact that during this phase the accretion rate does not reach the 
floor as it did in the previous intermediate stage and the oscilla- 
tions have progressively smaller amplitude while the accretion rate 
grows exponentially. 

It is worth pointing out that during each of these oscillations 
a considerable amount of matter falls towards the black hole and 
in a realistic scenario it is reasonable to expect that before this 
matter reaches the event horizon it will loose part of its poten- 
tial binding energy by increasing its temperature and by emitting 
electromagnetic radiation. In view of this, it is plausible to expect 
that the quasi-periodic accretion measured during our simulations 
could also be observed in the form of a quasi-periodic X-ray lu- 
minosity, as it is indeed the case in the quasi-periodic oscillations 
(QPO's) observed in the X-ray lumino sity of Low Mass X-ray Bi- 
naries (LMXB's) (van der Klis 200C). While the connection be- 



tween the two effects is very attractive, it should be remarked that 
the discs which are believed to be behind the quasi-periodic X- 
ray luminosity in LMXB's have much smaller rest-mass densities. 
Therefore, more detailed calculations need to be done before any 
firm conclusion on the connection between the two phenomenolo- 
gies can be drawn. 

All of the calculations shown in Figure [j] terminate when the 
accretion rate has reached a maximum value of ~ 3 x 10 3 M Q /s 
and the rest-mass of the torus has become only a few percent of the 
initial one (cf. the solid line of Figure ^|). During these very final 
stages of the instability the calculations become very difficult be- 
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cause of the exponential changes in the hydrodynamical quantities 
and Courant factors as small as 0.01 are needed to prevent the code 
from crashing. Soon after the accretion rate has reached its maxi- 
mum, it drops rapidly to very small values as a result of the almost 
complete disappearance of the torus (this final part of the evolution 
is not reported in Figure [j]). 

As mentioned above, the growth-rate of the runaway instabil- 
ity depends on the efficiency of the mass accretion process and on 
reaching a certain fractional change in the black hole mass. How 
rapidly this change takes place depends both on the strength of the 
initial velocity perturbation (as shown in Figure |l|) but also on the 
density of the accreted matter. To confirm this, we have performed 
simulations for the same initial perturbation but with different mass 
ratios Md/M. The results of these simulations are summarised in 
Figure ^ which shows the behaviour of the rest-mass accretion rate 
for models (b) and (c) in Table |l|. These models, we recall, have the 
same properties of model (a) but have been constructed using larger 
values of the polytropic constant k [cf. Eq. <Jsj>] - As expected, lower 
density toroidal neutron stars have proportionally smaller accretion 
rates (note that the floors in the mass flux m are different) and need 
longer timescales for the onset and development of the instability. 

Two interesting aspects of Figure Q need to be pointed out. 
Firstly, the timescale for the instability to set in has an almost lin- 
ear dependence on the mass ratio M t /M. This is an important de- 
tail as it reveals that what could be considered a "realistic" toroidal 
neutron star on the basis of the numerical simula t ions p erformed 



so far (Shibata & Uryu 



2000 



Ruffert & Janka 



20011, i.e. one 



with M t /M ~ 0.1 and with a level of perturbations of the order 
77 ~ 0.04, has a lifetime of roughly 0.2 s, if unstable to the runaway 
instability, (cf. also with Figure |l|). Secondly, when the mass accre- 
tion rate is generically low, the amount of matter accreted can be 
very small even over several tens of dynamical timescales. When 
this is the case, and as mentioned in Section the spacetime 
can be held fixed and the numerical calculations simplified. More 
importantly, this regime provides the possibility of distinguishing 
the dynamical response of the toroidal neutron star to perturbations 
from the development of the instability. We will exploit this possi- 
bility in the following Section which focuses on the investigation 
of the oscillation properties of toroiodal neutron stars. 



5.2 Fixed spacetime: quasi-periodic oscillations 

The dynamical response of the toroidal neutron star to perturba- 
tions provides information about the basic properties of this object 
in a strong gravitational field. While the details of these properties 
depend on the details of the gravitational field the torus is experi- 
encing, we expect some features to be generic and therefore to be 
present also in those circumstances in which the runaway instabil- 
ity is suppressed. To prevent the runaway instability from hiding the 
quasi-periodic response of the torus to perturbations, we have sim- 
ply suppressed the instability by maintaining the spacetime fixed. 

The most apparent consequence of this choice is that the 
torus is no longer dramatically accreted onto the black hole but re- 
mains at almost constant rest-mass for arbitrarily long times. This 
is shown in Figure |§| which displays the evolution of the rest-mass 
of the torus (normalized to its initial value) for both a dynamical 
(dashed line) and a fixed background spacetime, respectively. The 
solid line, in particular, refers to the model shown in Figure [l] with 
the same line style. Note that when the instability is fully devel- 
oped the toroidal neutron star has almost completely disappeared 
into the black hole, the remaining mass being just ~ 4% of the ini- 
tial one. Note also that the rest-mass evolution in a fixed spacetime 
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Figure 4. Rest-mass accretion rate in a fixed spacetime evolution. The solid, 
dotted and dashed lines correspond to r] = 0.06, 0.04 and 0.02, respec- 
tively. The inset shows the case r\ = 0.06 during the entire simulation (100 
orbital periods). 



is not exactly constant but shows a secular decay as a result of the 
small amounts of matter that are quasi-periodically accreted onto 
the black hole (cf. Figure ^). 

Figure ^ is the equivalent of Figure [l], showing the quasi- 
periodic accretion rate during the first 30 orbital periods for a sim- 
ulation in which the spacetime is held fixed. It is apparent that both 
evolutions have a qualitatively similar behaviour: after the toroidal 
neutron star has relaxed at t/t or b — 2, it starts accreting matter 
onto the black hole at quasi-periodic intervals which do not depend 
on the strength of the perturbation. The amplitude of the mass ac- 
cretion rate, on the other hand, does depend on the value of 77, pro- 
ducing larger amounts of accreted matter with increasingly larger 
values of 77. (Note that when r\ — 0.02, the mass accretion rate 
seems to drop to an almost constant value for 6 < t/t or b < 15- 
This is just the consequence of the logarithmic scale used and, as 
shown in Figure^, it is indeed possible to observe a periodicity also 
during this time interval.) 

Since the curves in Figure ^ are the result of numerical calcu- 
lations in which each period of oscillation requires several tens of 
thousands timesteps, the ability to reproduce this periodicity with 
such high precision is the result of the use of the accurate HRSC 
methods. In addition to this, the periodic behaviour does not seem 
to be altered even when observed over 100 orbital timescales as 
shown in the inset of Figure although some secular features ap- 
pear. Most notably, the mass accretion rate oscillates around val- 
ues that are increasingly larger. This is due to the fact that as the 
accretion proceeds, matter of increasingly larger rest-mass density 
reaches the cusp (the low-density, outer regions of the toroidal neu- 
tron star have already been accreted) and this therefore produces a 
small secular growth in the amplitude of the mass flux. 

The mass accretion rate is not the only quantity showing a 
periodic behaviour and indeed all of the fluid variables can be 
shown to oscillate periodically. In Figure ^ we show the time evolu- 
tion, over 30 orbital periods, of the central rest-mass density of the 
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Figure 5. Evolution of the central rest-mass density normalized to its ini- 
tial value. As in Figure |i| the solid, dotted and dashed lines correspond to 
r\ = 0.06, 0.04 and 0.02, respectively. The inset shows the case r\ = 0.06 
during the entire simulation (100 orbital periods). 



toroidal neutron star normalized to its initial value. The small in- 
set shows the evolution of the same quantity but for a much longer 
timescale and offers convincing evidence of the regularity of the 
oscillation. The inset should be contrasted with the evolution of the 
rest-mass density in a simulation with a dynamical spacetime (not 
shown here). In that case, in fact, the oscillations in the rest-mass 
density do not remain (roughly) constant but grow exponentially as 
the runaway instability develops. 

The periodic behaviour shown in Figure ^ has a simple inter- 
pretation: as a result of the initial perturbation, the toroidal neutron 
star acquires a linear momentum in the radial direction pushing it 
towards the black hole. When this happens, the pressure gradients 
become stronger to counteract the steeper gravitational potential 
experienced as the torus moves inward, thus increasing the central 
density and eventually pushing the torus back to its original po- 
sition. This is illustrated in more detail in the different panels of 
Figure ^, which show the rest-mass distribution at various times 
during one oscillation. 

More precisely, the sequence in Figure ^ shows that once the 
unperturbed toroidal neutron star [whose initial rest-mass density 
distribution is shown in the panel (a) of Figure ^] is subject to a 
radial velocity perturbation, it will start moving towards the black 
hole [panel (b)]. The existence of a potential barrier at r ~ 4.5, 
however, causes a compression of the matter that is approaching 
the black hole [panel (c)], giving rise to the first peak in the central 
density visible in Figure at t/t OI -t, ~ 0.45. Because the initial 
configuration is just marginally stable, a small fraction of the mat- 
ter in the torus will be pushed over the maximum of the potential 
barrier and generate the first maximum in the rest-mass accretion 
rate reported in Figure ^. The presence of a net mass flux onto the 
black hole can directly be appreciated through iso-density contours 
shown in the (r sin 6, r cos 6) planes of Figure^. Most notably, the 
lower density contours of panel (c) are closed on the event horizon 
and indicate therefore the presence of a thin channel of accreting 



matter that is linking the toroidal neutron with the black hole. (Be- 
cause of this correlation between the rest-mass density and the ac- 
cretion rate, the peaks in Figure^ are slightly advanced in time with 
respect to the corresponding peaks in Figure^.) As the compression 
increases, the pressure gradients become sufficiently strong to pro- 
duce a restoring force on the toroidal neutron star which is then 
pushed back, away from the black hole. The restoring effect is so 
efficient that the torus overshoots the original position [panel (a)] 
and moves outwards to larger radii [panel (d)]. When this happens, 
the central density decreases and the mass accretion rate drops to 
its floor value; both of these effects are reflected in the first minima 
of Figures ^ and ^. 

The dynamics of this process can also be followed by moni- 
toring the total energy of a fluid element at the edge of the torus, 
(ut)in. If, at a given time, this quantity becomes larger than the 
potential barrier at the cusp, W C us P , the corresponding fluid ele- 
ment will have sufficient kinetic energy to overcome the barrier but 
not sufficient angular momentum to sustain an orbital motion at the 
smaller radius at which it has been displaced. As a result, it will be 
forced to fall into the black hole, producing, after its free-fall time, 
a peak in the mass accretion rate. 

Once triggered, the behaviour described above will repeat it- 
self with great regularity (cf. the small insets of Figures ^ and 
and minimal numerical dissipation (which can be appreciated 
through the very small decay in the amplitude over 100 dynami- 
cal timescales) up until the numerical simulation is stopped or the 
toroidal neutron star has been entirely accreted by the black hole. 
During these oscillations, the pressure gradients act as a restoring 
force during the periodic transformation of the excess kinetic en- 
ergy (transferred with the initial velocity perturbation) into poten- 
tial gravitational energy and viceversa. As we shall comment on 
later, this is a first clue about the nature of these pulsations. 

It is worth noting that the introduction of a perturbation in 
the radial velocity will no longer guarantee that the specific angu- 
lar momentum is conserved during the time evolution. This can be 
most easily seen if we rewrite th e ^-component of the Euler equa- 
tions in the form [cf . Eq. (24) of ^Tawley et al. ( 1984a )] 



1 



d t t- 




(16) 



For an unperturbed toroidal neutron star, if — v e — = d^, 
and the right-hand-side of Eq. ( Jlrj ) is therefore zero, ensuring the 
conservation of the specific angular momentum. Clearly, this is no 
longer true when a small negative radial velocity is introduced and 
although small, this effect further favours the tiny spill of matter 
through the edge of the torus. 

A final comment should be made about initial models consist- 
ing of initially stable toroidal neutron stars ( AW < 0) and that are 
therefore fully contained in barotropic surfaces smaller than their 
Roche lobe. Also for these models we have performed a number 
of simulations investigating their behaviour for different initial per- 
turbations as well as for different initial masses. Overall, the be- 
haviour observed with these initial data is qualitatively similar to 
the one already discussed for marginally stable tori, i.e. also these 
models develop the runaway instability or show a quasi-periodic 
behaviour depending on whether the spacetime evolution is taken 
into account or not. The most significant difference observed is 
that smaller rest-mass accretion rates are generally produced for 
the same initial perturbation. This is simply due to the fact that in 
stable models a larger potential barrier is present at the cusp and is 
therefore increasingly more difficult for a fluid element to reach the 
black hole as a result of the initial perturbation. 
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Figure 6. Rest-mass density distributions of a perturbed toroidal neutron star at different times. The sequence illustrates the periodic behaviour during the first 
cycle of an oscillation excited with r) = 0.06. The iso-density contours plotted on the (r sin 0, r cos 9) planes of the different panels can be used to trace the 
motion of the matter and are the same for all of the panels (in particular they refer to p = 5.0xl0 12 ,l.x 10 13 , 2.5 X 10 13 , 5.0 X 10 13 , 7.5 X 10 13 , 1.0 X 10 14 , 
and 1.11 X 10 14 g/cm 3 , respectively). The times are expressed in terms of the orbital timescale, the spatial coordinates are in terms of the gravitational radius, 
and the rest-mass density in g/cm 3 . See main text for a comment on the sequence. 



5.2.1 Fourier analysis 

We have so far commented on the "quasi-periodic" behaviour of 
the hydrodynamical variables in response to the initial perturbation 
but we have not yet discussed how periodic is "quasi-periodic". 
We have therefore calculated the Fourier transforms of the relevant 
fluid quantities for a number of different models. As a good rep- 
resentative case, we plot in Figure |^ the power spectrum of the 
rest-mass accretion for models (a), (b) and (c) of Table |l|. The 
Fourier transform has been calculated with data obtained with an 
r\ — 0.06 perturbation and computed over a time interval going up 



to i/torb — 100. Note that larger values of r\ produce correspond- 
ingly higher peaks in the power spectra, but the data in Figure ^ has 
been suitably normalized to match the power in the first peak. 

There are two important features of Figure ^ that need to be 
pointed out. The first one is that all of the three power spectra shown 
consist of a fundamental frequency fo (228 Hz for the models con- 
sidered in Figure and a series of overtones (at 456 and 684 Hz, 
respectively) in a ratio which can be determined to be 1 : 2 : 3 : . . . 
and to an accuracy of a few percent. (Note that overtones higher 
than the second one have also been measured, although with much 
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Figure 7. Power spectrum of the mass accretion rate for models (a), (b), and 
(c). The values on the vertical axis have been suitably normalized to match 
the power in the fundamental mode and arbitrary units have been used. 

lower power than the one found in the first three peaks.) Indeed, the 
presence of at least three peaks can be detected also in the power 
spectra of basically all of the fluid variables as well as in the over- 
all displacement of the toroidal neutron star during the oscillations. 
Interestingly, the power spectra of some fluid variables, such as that 
of the 1/2 norm of the rest-mass density, seem to show peaks with 
much lower power also at intermediate frequencies, in particular 
at frequencies which are in a ratio 1:3/2:2:5/2:... with 
respect to the fundamental one. Because the energy of these peaks 
is very close to the background noise, it is not clear whether these 
modes correspond to physical modes or are due to numerical errors. 
The second important feature is that the peaks in the power spectra 
of the three models in Figure ^ have all the same frequencies, with 
differences below 0.1%. All of these properties are clearly suggest- 
ing that the quasi-periodic response observed is the consequence of 
some fundamental mode of oscillation of toroidal neutron stars and 
that, as for isolated neutron stars, it is probably independent even 
of the presence of a central black hole. 

Insight on the nature of these modes can be gained if one con- 
siders that all of the models reported in Figure M refer to toroidal 
neutron stars with fixed spatial dimensions and specific angular 
momentum, but with varying mass (cf. Table [l]). This has basi- 
cally been obtained by suitably rescaling the polytropic constant 
in the EOS. Another property shared by these models is that they 
all have the same initial average sound speed. We recall that for the 
ideal fluid configurations considered here, the local sound speed 
can be calculated as c s — 7(7 — l)p/[(y — l)p + jp\, which is 
effectively a constant for models with an initial density distribution 
given by Eq. (Q). As a result, it is not surprising that the peaks coin- 
cide in all of the models if the oscillations discussed so far should 
be associated to the p-modes (or acoustic modes) of the toroidal 
neutron star. However, proving that the periodic oscillations ob- 
served in our simulations can be interpreted as p-modes and vali- 
dating that the presence of modes at intermediate frequencies is not 
a numerical artifact requires a detailed perturbative analysis and is 
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Figure 8. Dependence of the fundamental frequency fo on the average den- 
sity in a sequence of models having the same mass ratio, M^/M = 0.1. 
The data refers to simulations with initial perturbations of strength 77 = 
0.06. 



beyond the scope of this paper. Work is now in progress and pre- 
liminary results in this direction seem to confirm the hypothesis 
that these oscillations represent the vibrational modes of relativis- 
tic toroidal neutron stars having time-va rying pressure gradients as 
the restoring force (Rezzolla et al. 2002). 



Since the peaks in the power spectra seem related to the av- 
erage sound speed, the dependence of the fundamental frequency 
on the properties of the toroidal neutron star needs to be investi- 
gated along sequences different from the one consider in Figure ^. 
As an example, we have reported in Figure [§] the fundamental fre- 
quency fo as a function of the average rest-mass density inside the 
torus. The data refer to the evolution of models (e)-(h) in Table [l], 
that have different dimensions, specific angular momenta and poly- 
tropic constants, but all have the the same mass ratio M t /M (cf. 
Table [l]). Figure ^ shows that an evident correlation exists between 
the fundamental frequency and the logarithm of the average rest- 
mass density. A fit to the data indicates that this correlation is (cf. 
straight line in Figure ^) 



fo = (122.55 log (p) - 1201.67) Hz . 



(17) 



where (p) is expressed in cgs units. This expression is important as 
it suggests that a systematic study of these oscillations for different 
initial models of toroidal neutron stars is possible. Furthermore, it 
represents a first step towards a relativistic disc-seismology analy- 
sis for massive and vertically extended tori in General Relativity, 
in anal ogy to the o ne extensively dev e loped for geometrically thin 
discs ( frCato| |200l|; jPerez et al.||l997t ISilbergleit et~al| |200l[ R o- 



driguez et al. |2002[ ) 



5.2.2 Linear and nonlinear regimes 

All of the quasi-periodic behaviour discussed so far is the conse- 
quence of the finite size perturbations that have been introduced in 
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Figure 9. Averaged maximum rest-mass density normalized to the central 
one in the toroidal neutron star as a function of the perturbation strength. 
The small inset shows a magnification of the behaviour in the lineal' regime 
and the solid line shows a very good linear fit to the data. The data refers to 
model (a) and has been averaged over 10 orbital timescales. 



the initial configuration. Within this approach a linear regime is ex- 
pected in which the response of the toroidal neutron star is linearly 
proportional to the perturbation introduced, and a nonlinear regime 
when this ceases to be true. The strength of the perturbation which 
marks the transition between the two regimes can be estimated from 
Figure ^, where we show the averaged maximum rest-mass density 
normalized to the central one in the toroidal neutron star. 

A rapid look at Figure ^[ in fact, reveals the presence of both 
the linear and nonlinear regimes with the first one being shown 
magnified in the inset, and where the solid line shows the very good 
linear fit to the data. The transition between the two regimes seems 
to occur for r\ ~ 0.05, with the nonlinear regime producing maxi- 
mum amplitudes that are ~ 35% larger than the initial one. A care- 
ful analysis of the behaviour of the fluid variables shows that for 
perturbations with strength r\ > 0.05, some of the kinetic energy 
which is confined to the lower order modes in the linear regime, 
tends to be transferred also to higher order modes (As remarked 
above, while in Figure^ we have reported only the first three peaks, 
the power spectra show the presence of higher order overtones, with 
up to the seventh one being clearly visible). The nonlinear coupling 
among different modes and the excitation of higher order overtones 
is often encountered in Nature where it serves to redistribute the 
excess kinetic energy before the production of shocks. In practice, 
the nonlinear coupling deprives of energy the fundamental mode 
(which is the one basically represented in Figure^) and is therefore 
responsible for the decay of (p ma x/pc) for r\ > 0.1. Interestingly, 
when analysed in terms of the power spectra, this effect shows a 
very distinctive behaviour. As the nonlinear mode-mode coupling 
becomes effective, the amount of power in the fundamental mode 
becomes increasingly smaller as the strength of the perturbation 
is increased. At the same time, the conservation of energy trans- 
fers power to the overtones, with the first ones reaching amplitudes 



comparable to the fundamental one and with the high order ones 
becoming more and more distinct from the background. 

Determining the transition to the nonlinear regime is impor- 
tant to set an approximate upper limit on the amplitude of the oscil- 
lations and, as we will discuss in the following Section, it is relevant 
when estimating the emission of gravitational waves. It should also 
be noted that in the parameter range for r\ in which we have per- 
formed our calculations (i.e. r\ £ [0.001, 0.12]), the peak frequen- 
cies in the power spectra have not shown to depend on the values 
used for r\. This is of course consistent with them being fundamen- 
tal frequencies (and overtones). 

A final comment should be made on the minimum value of 
the perturbation parameter r\ which is sufficient to produce the 
quasi-periodic behaviour described in this Section. On the basis of 
the continuum equations one expects that this minimum value is 
strictly larger than zero. However, we have performed simulations 
of marginally stable configurations (AW = 0) with r\ = and 
observed much of the phenomenology described above even if at 
rather minute amplitudes and with a larger numerical noise. This 
result, which has been e ncountered also in other accurate simula- 
tions of relativistic stars (Font et al. |2002| ), is not entirely surprising 
and is simply indicating that in the absence of prescribed perturba- 
tions, even the small truncation error introduced in the construction 
of the initial configuration is sufficient to excite the pulsations in 
these modes of vibration. 



6 GRAVITATIONAL WAVE EMISSION 

Despite our analysis has been restricted to axisymmetry, two sim- 
ple considerations suggest that strong quasi-periodic gravitational 
waves should be expected together with the quasi-periodic accre- 
tion. The first of these considerations is that these oscillating tori 
undergo large and rapid variations of their mass quadrupole mo- 
ment which, we recall, can be calculated as 



, 3 2 
P I -z 



1 



r 4 drdz , 



(18) 



where z = cos 8. It is therefore reasonable to expect that gravita- 
tional waves, with a strain proportional to the second time deriva- 
tive of Eq. should be produced during such oscillations^ The 
second consideration is suggested by expression (|l8|), which shows 
that a configuration with toroidal topology and in which the rest- 
mass density has a maximum away from the origin will naturally 
have a large mass quadrupole simply because of the product pr 4 in 
the integrand of this equation. This should be contrasted with what 
happens for stars with spherical topology and which have instead 
the largest densities at the centre. Both of these arguments justify 
the intense gravitational radiation and the large signal-to-noise ra- 
tio we have calculated using the Newtonian quadrupole approxi- 
mation. Before discussing this in detail, however, it is useful to re- 
mind that since we are not solving the Einstein field equations, we 
are unable to account for that part of the gravitational radiation that 
is emitted by the black hole itself as a result of the quasi-periodic 
mass accretion. For the same reason we cannot estimate the amount 
of gravitational radiation that will be captured by the black hole and 
will not reach null infinity. Work is now in progress to calculate 
also this part of the radiative field using an approach in which the 

2 It is worth underlying that a deviation from axisymmetry is only a suf- 
ficient condit ion for the emission o f gravitational waves and not, as some- 



times stated ( |vlineshige et al. 2002 1, a necessary condition 
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fluid evolution is used as a "source" for a perturbative form of the 
Einstein equations. 

The presence of an azimuthal Killing vector has two important 
consequences. Firstly, the gravitational waves produced by these 
axisymmetric oscillating tori will carry away energy but not an- 
gular momentum, which is a conserved quantity in this spacetime. 
This is to be contrasted with the gravitational wave emission from 
non-axisymmetric perturbations in a massive torus orbiting a black 
ho le, whose strength has been estimated in recent papersj van Put- 



tcn 2001e ^r Vlineshige et al. 2002 ). Secondly, the gravitational 
waves produced will have a si ngle polarization s tate ( i.e. the "plus" 
one in our coordinate system, Kochanek et al. (199C)), so that the 



transverse traceless gravitational field is completely determined in 
terms of its only nonzero trans verse and traceless (TT) i ndependent 
components, h^J — —h^ (Zwerger & Mullet 1997). Adopting 
then Newtonian quadrupole approximation, we can calculate the 
gravitational waveform h TT (t) observed at a distance R from the 
source in terms of the quadrupole wave amplitude Afg (see also 
Zwerger & Miiller 1997) 



R) 



R 



(19) 



where F + — F + (R, 8, (j>) is the detector's beam pattern function 
and depends on the orientation of the source with respect to the 
observer. As customary in these calculations, we will assume it to 
be optimal, i.e. F+ = 1. The wave amplitude A^q in Eq. ( |T^ ) is 
simply the second time derivative of the mass quadrupole moment 
and can effectively be calculated without taking time derivatives, 
which are instead replaced by spatial derivatives after exploiting 
the continuity and the Eul er equations (Finn 



2001 



Rezzolla et al 



198? ; Blanchet et al 



A' 



dt 2 



1999|) to give 

r/ Q 2 
V r V (6z - 



1)+V0V (2 — 3z )—V c f,V 4 



-62y {y r v r )(vev e )(l — z 2 



9 * fl 2 



drdz . 



(20) 



where k = /\/l5, and <3> is the gravitational potential. 

Since Eq. ( |2c| ) is intrinsically Newtonian, it brings up two subtle 
issues when evaluated within a relativistic context. These basically 
have to do with the definition of the radial coordinate and with the 
definition of the gravitational potential appearing in Eq. (|2c|). We 
have here opted for a pragmatical approach and treated r as the 
Schwarzschild radial coordinate and computed the gravitational po- 
tential in terms of the radial metric function as $ = (1 — g rr )/2, 
which is correct to the first Post-Newtonian (PN) order. 

To validate the correct implementation of the integral in 
Eq. (j2C|), we have considered a stationary torus in stable equilib- 
rium, (i.e. A Win < 0) and without any perturbation besides the one 
introduced by the truncation error (i.e. rj — 0). Under these circum- 
stances, no gravitational radiation can be produced and the terms in 
the square brackets of Eq. ( pfj[ ) should therefore compose to give 
zero identically. When computed numerically, we have found that 
the sum of these terms is effectively very small and of the order 
~ 10 -2 . This small residual in the integrand is due to approxi- 
mations mentioned above (i.e. the use of the Schwarzschild radial 
coordinate and the first PN approximation to the gravitational po- 
tential) and should be interpreted, in practice, as a consequence of 
the fact that our tori are highly relativistic objects, whose equilib- 
rium is not exactly given by the balancing of the Newtonian terms 
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Figure 10. Time evolution of the second time derivative of the mass 
quadrupole, computed for model (a). The solid, dotted and dashed lines 
correspond to r\ = 0.06, 0.04, 0.02, respectively. The inset corresponds to 
the case rj = 0.06 and spans 100 orbital timescales. 



in the square brackets. Because of this residual, however, the wave 
amplitude Afo will not average to zero over an oscillation but will 
have a net offset. We account for this by removing the overall resid- 
ual in the evaluation of the wave amplitude (see also Dimmelmeier 
et al. 2002 for the discussion of an analogous technique). 

Figure |To| shows the time evolution of the wave amplitude Af 2 
computed for model (a) via the expression (|2o|). The different line 
types refer to 77 = 0.06, 0.04 and 0.02, respectively. The com- 
puted gravitational waveform exhibits the same periodic behaviour 
discussed in the previous Sections for several fluid variables and 
shows oscillations that are in phase with the ones observed in the 
rest-mass density (cf. Figure p|). Furthermore, and as one would 
expect, the wave amplitude scales with the strength of the initial 
perturbation. Note that in the case of a simulation with a dynami- 
cal spacetime (not shown in Figure |Io|), the gravitational waveform 
does not maintain a constant in time amplitude but, as the run- 
away instability develops, the variations in Af ( 2 become increas- 
ingly large, with an exponential growth rate that matches the one 
observed in the density evolution. 

Since the mass quadrupole and its second time derivative are 
linear in the rest-mass density [cf. Eq. @] and the latter exhibits 
a quasi-periodic behaviour upon perturbations, it is natural to ex- 
pect the same linear dependence to be present also in terms of the 
mass ratio Mt /M. To verify this we have computed Afo with an 
initial perturbation rj = 0.06 for models (a), (b), and (c) that, we 
recall, differ only for their mass Mt. The results of this analysis 
are reported in Fig. which shows the behaviour of Afo over one 
representative period of oscillation and for the four models. As ex- 
pected, the scaling of the amplitude is linear with Mt/M and this 
is more clearly shown in the inset where the wave amplitude [in- 
cluding the data for model (d)] is fitted linearly with Mt/M. 

Using Eq. <[l^), it is possible to derive a phenomenological 
expression for the gravitational waveform that could be expected 
as a result of the oscillations induced in the toroidal neutron star. 
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Figure 11. Dependence of the wave amplitude Alj? on the mass of the 
torus. The solid, dotted and dashed lines correspond to models (a), (b) and 
(c), respectively. The data refers to an initial perturbation of strength rj = 
0.06 and the inset shows a linear fit (solid straight line) of the data, including 
the one for model (d). 



Restricting our analysis to the linear regime for which a simpler 
scaling is possible, we can express the transverse traceless gravi- 
tational wave amplitude for a source in the Galaxy in terms of the 
relevant parameters in our problem 



h TT ~ 2.2 x 10" 



___ 
0.04 



Mt 



0.1 M 2 ._ 



10 Kpc 



(21) 



where we have defined M2.5 = M/(2.5 Mq). 

Overall, expression ( |_"l| ) shows that, already in the linear 
regime for the perturbations, a non-negligible gravitational wave 
amplitude can be produced by an oscillating toroidal neutron star 
orbiting around a black hole. This amplitude is indeed compara- 
ble with the average gravitational wave amplitu de computed in the 



case of core collapse in a super nova explosion (Zwerger & Mullet 



1997 ; Dimmelmeier et al 



2002 ) and can become stronger for larger 
perturbations or masses in the torus. 

Of course, the large wave amplitudes suggested by expression 
( p"jj ) are referred to a galactic source and would become three orders 
of magnitude smaller for a source located at the edge of the Virgo 
cluster (i.e. at about 20 Mpc). What is important to bear in mind 
from expression ( ^l[ ) is that oscillating toroidal neutron stars can 
be sources of gravitational waves as strong or stronger than a core 
collapse and that could occur with comparable event rates^J. This 
notion then serves as a useful normalization for estimating their 
relevance. It should also be noted, however, that a strong gravita- 
tional wave signal is just a necessary condition for the detectability 
of the gravitational wave emission from toroidal neutron stars and 
that additional conditions, such as a sufficiently high event rate or 

3 Of course not all of the collapsing iron cores will produce a black hole 
sourrou nded by a massive torus. On the other hand, and as discussed in 
Section there is a multiplicity of scenarios in which toroidal neutron 
stars could be produced. 



a good matching with the sensitivies of the detectors, need to be 
met. In the following Section we will discuss these issues in more 
detail and evaluate the detectability of these potential new sources 
of gravitational waves. 

6.1 Detectability 

To assess the detectability of toroidal neutron stars as sources of 
gravitational radiation we have computed the characteristic gravi- 
tational wave frequency and amplitude, as well as the correspond- 
ing signal-to-noise ratio for the interferometric detectors that will 
soon be operative. More specifically we have first computed the 
gravitational waveform in the frequency domain as the Fourier 
transform of the traceless transverse waveform in the time domain 



Hf) 



2-Kift . TT 



(t) dt , 



(22) 



where h TT (t) is calculated according to Eq. jl9| ) and where we 
have considered the gravitational wave amplitude computed over a 
fixed spacetime only. (Hereafter we will indicate h TT (t) simply as 
h(t).). While in principle the integral in Eq. ( ^_| ) is over an infinite 
time interval, in practice h(t) is nonzero only over a finite inter- 
val Tijf c . If the initial model chosen for the toroidal neutron star is 
not a stable one, this time interval is simply set by the timescale 
over which the runaway instability develops. If, on the other hand, 
the initial model is stable, the timescale over which a gravitational 
wave signal is produced can be considerably longer and is basically 
set by the time over which the star is able to survive, for instance, 
against non-axisymmetric instabilities. Hereafter, as a representa- 
tive timescale for our "realistic" toroidal neutron star model (cf. 
Section 5.1) we will assume Tus e ~ 0.2 s, while longer/shorter 



timescales will be assumed for models with smaller/larger initial 
perturbations (cf. Table ^|). 

Given a detector whose response has a power spectral density 
Sh(f), it is then use ful to calculate characteristic frequency of the 
signal ([Thorne||l987|) 



(\Kf)\ 2 ) 
s h (f) 



fdf 



(Hf)\ 2 ) 

Sh(f) 



df 



(23) 



where (\h(f)\ 2 ) denotes an average over randomly distributed an- 
gles, that we have simply approximated as {\h(f)\ 2 ) ~ \h(f)\ 2 . 
The characteristic frequency provides a representative measure of 
where, in frequency, most of the signal is concentrated and is there- 
fore relevant when the gravitational wave signal has a rather broad 
spectrum in frequency, as is the case in a gravitational core col- 
lapse. In the case of an oscillating toroidal neutron star, on the other 
hand, the signal is basically emitted at the fundamental frequency 
of oscillation (cf. Table ^), increasing the detectability. 

Once the characteristic frequency is known, it can be used to 
determine the characteristic amplitude as 



Sh(fc) 

Sh(f) 



WW) fdf 



1/2 



(24) 



It worth remarking that when the characteristic frequency fits well 
in the minima of the sensitivity curves, the weight Sh(fc) / Sh(f) 
in the integral of ( p4| ) can significantly increase h c when compared 
to h TT . This is indeed what happens for the toroidal neutron stars 
considered here (cf. Table Q). 

A direct comparison of characteristic amplitude with the root- 
mean-square strain noise of the detector 



firms = v fSh(f) 



(25) 
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Table 2. Characteristic properties for the detection of the gravitational wave signal emitted by a toroidal neutron star. The first two columns show the mass of 
toroidal neutron star normalized to the black hole one, as well as the strength of the initial perturbation r). The following nine columns report the characteristic 
frequency, the characteristic wave strain and the signal-to-noise ratio computed for the three detectors LIGO I, LIGO II and EURO, assuming a galactic 
distance for a source detected by LIGO I, and an extragalactic distance for a source detected by LIGO II and EURO. The last column reports the characteristic 
lifetime for the existence of toroidal neutron stars that are unstable to the runaway instability. 
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noise of LIGO I and LIGO II, but the strain curves of VIRGO (that 
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Figure 12. Characteristic wave amplitudes for a perturbed toroidal neutron 
star with rj = 0.04. These amplitudes have been computed using the strain 
noise estimated for LIGO I for a source located at a distance of 10 Kpc and 
for LIGO II for a source at a distance of 20 Mpc. The numbers can also be 
effectively compared with the sensitivity curves of VIRGO and EURO that 
are similar in this frequency range. Different points refer to different mass 
ratios, with triangles indicating Mt/M = 0.25 and squares Mt/M = 0.1. 



finally determines the signal-to-noise ratio at the characteristic fre 
quency as 

S_ h, 
N 



f \ ■ (26) 

In Figure [l2| we show the characteristic wave amplitude for 
sources located at a distance of R — 10 Kpc and R = 20 Mpc, 
as computed for two different values of the toroidal neutron star 
mass. These amplitudes have been computed for the expected strain 



is simil ar in this freque ncy range, Damour et al (2001)) and of 
EURO (Winkler 2001) have also been reported for comparison. 
Interestingly, with small initial perturbations (77 = 0.04) and mass 
ratios (Mt/M = 0.1), the computed characteristic amplitudes can 
be above the sensitivity curves of LIGO I for sources within 10 
Kpc and above the sensitivity curve of EURO for sources within 20 
Mpc. Both results suggest that a toroidal neutron stars oscillating 
in the Virgo cluster could be detectable by the present and planned 
interferometric detectors. 

Summarized in Table |^ are the basic properties of the grav- 
itational wave signal emitted by a toroidal neutron star. Most no- 
tably, we report: the characteristic frequency, the characteristic 
wave strain and the signal-to-noise ratios as computed for LIGO I, 
LIGO II and EURO, as well as the timescale over which the signal 
has been computed. All of these quantities refer to models with dif- 
ferent initial perturbations and located at either 10 Kpc (for LIGO 
I) or 20 Mpc (for LIGO II and EURO). 

The values for the signal-to-noise ratios reported in Table ^ are 
already interesting, but could become larger in at least three differ- 
ent ways. Firstly, and as remarked above, the signal strength com- 
puted does not include the exponential growth in the gravitational 
waveform that would accompany the runaway instability and that 
would provide an important contribution to the overall characteris- 
tic amplitude, as it is the case in binary mergers. As an illustrative 
example, consider that if the last 10 ms (i.e. roughly the last 5 os- 
cillations) before the disappearence of an unstable torus at 10 Kpc 
with Mt/M = 0.1 and an initial perturbation with rj = 0.02 were 
taken into account, they would yield a final S/N ~ 90. This is to 
be contrasted with the corresponding S/N = 11.9, obtained when 
the oscillations are considered on a fixed spacetime (cf. Table |^). 
Secondly, we here have assumed the lifetime of the tori to be lim- 
ited by the runaway instability to avoid the uncertainties related to 
the very existence of the instability when the specific angular mo- 
mentum is not constant. On the other hand, a toroidal neutron star 
oscillating for a timescale longer than the one assumed in Table ^ 
will have a proportionally stronger signal even when the exponen- 
tially growing phase is neglected. Again, as an example, it is useful 
to consider that the model yielding S/N = 11.9 over 0.2 s would 
produce a S/N — 16.3 if ruf e — 0.4 s. Thirdly, an oscillating 
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toroidal neutron star stable to the runaway instability would proba- 
bly be subject to viscous or magnetic driven non-axisymmetric in- 
stabilities on timescales longer than the ones discussed here. Once 
the non-axisymmetric deformations are fully developed the torus 
would then have a mass quadrupole with much larger time varia- 
tions, losing amounts of energy (and angular momentum) in gravi- 
tational waves larger than the ones computed here. 

The final point to be addressed is the rate at which the emis- 
sion of gravitational waves from toroidal neutron stars could be 
detected. Although a realistic estimate of this rate is very difficult 
since very little is still known about the formation of toroidal neu- 
tron stars, it is reasonable to expect that these objects will be pro- 
duced in a significant fraction of the events leading to core collapse 
in supernova explosions, binary neutron star mergers and tidal dis- 
ruption of neutron stars orbiting a black hole. Because in a volume 
comprising the Virgo cluster all of these sources are among the 
most promising ones, overall our results indicate that toroidal neu- 
tron stars could potentially be new sources of gravitational radiation 
and certainly suggest a more accurate analysis. 



7 CONCLUSIONS 

We have performed general relativistic hydrodynamics simulations 
of axisymmetric massive and compact tori orbiting a Schwarzschild 
black hole with a constant specific angular momentum. These ob- 
jects have hydrostatic and hydrodynamical properties (i.e. large 
masses in small volumes with central rest-mass densities reaching 
almost nuclear matter density) that make them behave effectively 
as neutron stars, while possessing a toroidal topology. Our attention 
has been concentrated on two different aspects of the dynamics of 
these tori. 

The first one was focussed on the suggestion that the toroidal 
neutron stars might be dynamically unstable to the runaway insta- 
bility only if suitably chosen initial data were prescribed. To inves- 
tigate this we have used a time-dependent numerical code that inte- 
grates the general relativistic hydrodynamics equations on a curved 
background using HRSC methods. The evolution of the spacetime 
is an essential feature for the development of the instability and 
has been modelled through a sequence of stationary Schwarzschild 
spacetimes differing only in the total mass content, which is com- 
puted in terms of the total rest-mass accreted onto the black hole 
at each instant. The conclusion reached is that, at least for constant 
specific angular momentum tori whose self-gravity is neglected, the 
initial Roche lobe overflow is not a necessary condition for the de- 
velopment of the instability, which represents a natural feature of 
the dynamics of these objects. Very recent numerical work shows 
that this conclusion does not nece ssarily hold when the s pecific an- 
gular momentum is not constant (Font & Daigne 2002b). 



The second aspect was focussed on the dynamical response 
of these relativistic tori to the perturbations that are expected to be 
present after the catastrophic events that lead to their formation. 
Upon the introduction of suitably parametrized perturbations, the 
toroidal neutron stars have shown a regular oscillatory behaviour 
resulting both in a quasi-periodic variation of the mass accretion 
rate as well as of the rest-mass distribution. This response has been 
interpreted in terms of the excitation of oscillations modes which 
could be associated with the p-modes of toroidal neutron stars. 
These modes, which have been detected both in their fundamen- 
tal frequencies as well as in their overtones, depend systematically 
on the average density of the tori, and a disco-seismologic analysis 



could provide important information on the physical properties of 
these toroidal neutron stars. 

High rest-mass densities together with a toroidal topology are 
the basic properties that yield large mass quadrupoles for these 
toroidal neutron stars. As a consequence of the excitation of os- 
cillations, the mass quadrupoles are induced to change rapidly and 
intense gravitational radiation is thus produced. Estimates made 
within the Newtonian quadrupole approximation have shown that 
strong gravitational waves can be produced during the short life- 
time of these tori. In particular, the gravitational radiation emitted 
by these source is comparable or larger than the one that is ex- 
pected during the gravitational collapse of a stellar iron core, with 
a rate of detectable events which could also be larger given the 
variety of physical scenarios leading to the formation of a massive 
torus orbiting a black hole. Overall, the strength of the gravitational 
waves emitted and their periodicity are such that signal-to-noise ra- 
tios ~ 0(1) — 0(10) can be reached for sources at 20000 or 10 
Kpc respectively, making these new sources of gravitational waves 
detectable and potentially important. 

The results reported here are a first step towards the under- 
standing of the dynamics of toroidal neutron stars and call for a 
number of natural extensions and improvements. In particular, we 
plan to extend the simulations to the more realistic scenario of non- 
constant angular momentum tori, in order to find out if our conclu- 
sions regarding both the oscillation properties of toroidal neutron 
stars and the large amplitude of the associated gravitational wave 
emission still hold. In addition, it is interesting to improve, using 
more accurate approaches, the calculation of the gravitational radi- 
ation from the oscillations as well as from the runaway instability 
in t he region of the paramet er space where the instability could ex- 
ist (Font & Daigne 2002b). Equally interesting is to include the 



effects introduced by self-gravity of the torus and determine the 
quantitative differences that will be encountered in this case. Last 
but not least, it is worth analyzing how the presence of a Kerr black 
hole would modify the present results, determining, in particular, 
the dependence of this phenomenology on the spin of the central 
black hole. 

The results obtained here also provide promising suggestions 
for further work in a number of research areas different from the 
ones considered here. Firstly, they indicate that more needs to be 
investigated about the oscillation properties of relativistic tori and 
geometrically thick discs. Secondly, they suggest the possibility 
that the phenomenology observed in the quasi-periodic X-ray lu- 
minosity existing in LMXB's might be related, in some form, to 
the qu asi-periodic accretion resulting from oscillation mo des of the 
discs (Markovic & Lamb||l998i |Psaltis & Normar] |2000j W agoner 
et al. |2002j ). Thirdly, and perhaps most importantly, these results 
show that new and unexpected sources of gravitational radiation 
could exist and might be observed when the new detectors of grav- 
itational radiation become fully operative in the near future. 
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